Opis konkursu

Zadanie polega na prognozie sprzedaży produktów z różnych dziedziń w sieci sklepów Favorita w Ekwadorze.

Zestaw treningowy i testowy zawierają następujące dane:

• sales - tylko w zestawie treningowym, ilość sprzedanych produktów, ułamki pochodzą z produktów na wagę
• store_nbr - numer sklepu, w sumie jest 54 różnych sklepów
• family - dziedzina produktów np. pieczywo, nabiał, książki, etc., w sumie jest 33 różnych dziedzin
• onpromotion - ilość produktów na promocji
• date - data korespondująca z wyżej wymienionymi danymi danego dnia
• id - indeks

Ponadto do wykorzystania są następujące zestawy danych:

• stores - informacje o lokalizacji sklepów - miasto, rewir/prowincja, klasa sklepu - te informacje nie zostały wykorzystane w modelowaniu
• oil - cena ropy danego dnia, od której gospodarka Ekwadoru jest silnie zależna
• holiday_events - informacje dotyczące świąt i wydarzeń takie jak czy święto jest obchodzone w całym kraju, czy tylko w pewnych rejonach, czy święto/wydarzenie zostało przeniesione na inny dzień

Zestaw testowy zawiera informacje na 15 dni dla których należy przewidzieć ilości sprzedaży produktów każdej dziedziny w każdym ze sklepów.

Opis konkursu zawierał dwie dodatkowe informacje:

• wypłaty odbywają się w połowie miesiąca
• w połowie kwietnia 2016 roku zaszło trzęsienie ziemi) w północno zachodniej części Ekwadoru, którego efekty mogą być widoczne w sprzedaży w niektórych sklepach (zniszczona infrastruktura, pomoc charytatywna, priorytetyzowanie niezbędnych produktów)

Link do konkursu na Kaggle i oryginalnego opisu.

Spis treści:

• Eksploracja danych
• Analiza czasowa
• Przygotowanie zestawów treningowych
• Modelowanie
• Wnioski

import pandas as pd
import numpy as np

import math
import matplotlib.pyplot as plt
import seaborn as sns
import textwrap
import json

import plotly.express as px
from plotly.subplots import make_subplots
import plotly.graph_objs as go
import plotly.io as pio
from IPython.display import Image

from sklearn.linear_model import LinearRegression
from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess

from datetime import datetime, timedelta
import warnings

# Prosta funkcja do wyświetlania wykresów plotly jako obraz rastrowy:
PNG_PLOTS = True

def check_png(fig, rule=PNG_PLOTS):
    if rule:
        img_bytes = fig.to_image(format="png")
        fig = Image(img_bytes)
        return fig
    else:
        return fig

df = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')
holidays = pd.read_csv('holidays_events.csv')
oil = pd.read_csv('oil.csv')
stores = pd.read_csv('stores.csv')
transactions = pd.read_csv('transactions.csv')

df

	id	date	store_nbr	family	sales	onpromotion
0	0	2013-01-01	1	AUTOMOTIVE	0.000	0
1	1	2013-01-01	1	BABY CARE	0.000	0
2	2	2013-01-01	1	BEAUTY	0.000	0
3	3	2013-01-01	1	BEVERAGES	0.000	0
4	4	2013-01-01	1	BOOKS	0.000	0
...	...	...	...	...	...	...
3000883	3000883	2017-08-15	9	POULTRY	438.133	0
3000884	3000884	2017-08-15	9	PREPARED FOODS	154.553	1
3000885	3000885	2017-08-15	9	PRODUCE	2419.729	148
3000886	3000886	2017-08-15	9	SCHOOL AND OFFICE SUPPLIES	121.000	8
3000887	3000887	2017-08-15	9	SEAFOOD	16.000	0

3000888 rows × 6 columns

Eksploracja danych

# funkcja do wyświetlania informacji o ramkach danych

def df_overview(dataframe, modulo=5, head_rows=3, zip_nans=True):
  shape_rows = dataframe.shape[0] # ile wierszy
  shape_cols = dataframe.shape[1] # ile wszystkich kolumn
  shape_cols_rest =  shape_cols % modulo # ile kolumn po dzieleniu (ogon)

  main = dataframe.iloc[:,:-shape_cols_rest] # df z głównymi kolumnami
  rest = dataframe.iloc[:,-shape_cols_rest:] # df z ogonem kolumn

  main_cols = np.array(main.columns.tolist()).reshape(int(main.shape[1] / modulo), modulo) 
  # nazwy głównych kolumn w grupach
  rest_cols = rest.columns # nazwy kolumn ogon

  for i in main_cols: # do wyświetlenia kolumn z ilorazu
    head = dataframe[i].head(head_rows) # góra outputu

    uniques = [dataframe[j].nunique() for j in head.columns] # ile unikatowych wartości

    nans =  [dataframe[j].isna().sum() for j in head.columns] # ile NaN
    nans_perc = [round(dataframe[j].isna().sum() / shape_rows  * 100, 1) for j in head.columns ] # % NaN

    nans_and_perc = [(nans[i], str(nans_perc[i])+' %') for i in range(len(nans))] 
    # kombinacja obu żeby była jeden wiersz

    types = [dataframe[dataframe[j].notnull()][j].dtype for j in head.columns] # jakie typy


    if zip_nans: # czy chcemy skrócone NaN
      head.loc[len(head)] = types
      head.loc[len(head)] = uniques
      head.loc[len(head)] = nans_and_perc

      head = head.rename(index={int(len(head)-3): 'Data type',
                        int(len(head)-2): 'Unique values',
                        int(len(head)-1): 'NaNs',
                        })

    else: # czy chcemy pełne NaN
      head.loc[len(head)] = types
      head.loc[len(head)] = uniques
      head.loc[len(head)] = nans
      head.loc[len(head)] = nans_perc

      head = head.rename(index={int(len(head)-4): 'Data type',
                        int(len(head)-3): 'Unique values',
                        int(len(head)-2): 'NaN total',
                        int(len(head)-1): 'NaN %'
                        })

    display(head) # wyświetla kolumny z ilorazu


  for i in [rest_cols]: # do wyświetlenia kolumn z reszty, robimy [list] żeby działało tak samo
    head = dataframe[i].head(head_rows) # góra outputu

    uniques = [dataframe[j].nunique() for j in head.columns] # ile unikatowych wartości

    nans =  [dataframe[j].isna().sum() for j in head.columns] # ile NaN
    nans_perc = [round(dataframe[j].isna().sum() / shape_rows  * 100, 1) for j in head.columns ] # % NaN

    nans_and_perc = [(nans[i], str(nans_perc[i])+' %') for i in range(len(nans))] 
    # kombinacja obu żeby była jeden wiersz

    types = [dataframe[dataframe[j].notnull()][j].dtype for j in head.columns] # jakie typy
    if zip_nans: # czy chcemy skrócone NaN
      head.loc[len(head)] = types
      head.loc[len(head)] = uniques
      head.loc[len(head)] = nans_and_perc
      head = head.rename(index={int(len(head)-3): 'Data type',
                        int(len(head)-2): 'Unique values',
                        int(len(head)-1): 'NaNs',
                        })

    else: # czy chcemy pełne NaN
      head.loc[len(head)] = types
      head.loc[len(head)] = uniques
      head.loc[len(head)] = nans
      head.loc[len(head)] = nans_perc
      head = head.rename(index={int(len(head)-4): 'Data type',
                        int(len(head)-3): 'Unique values',
                        int(len(head)-2): 'NaN total',
                        int(len(head)-1): 'NaN %'
                        })
    print(dataframe.shape)
    display(head) # wyświetla kolumny z reszty

print('train:')
display(df_overview(df, 10))
print('test:')
display(df_overview(test, 10))
print('holidays:')
display(df_overview(holidays, 10))
print('oil:')
display(df_overview(oil, 10))
print('stores:')
display(df_overview(stores, 10))
print('transactions:')
display(df_overview(transactions, 10))

train:
(3000888, 6)

	id	date	store_nbr	family	sales	onpromotion
0	0	2013-01-01	1	AUTOMOTIVE	0.0	0
1	1	2013-01-01	1	BABY CARE	0.0	0
2	2	2013-01-01	1	BEAUTY	0.0	0
Data type	int64	object	int64	object	float64	int64
Unique values	3000888	1684	54	33	379610	362
NaNs	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)

None

test:
(28512, 5)

	id	date	store_nbr	family	onpromotion
0	3000888	2017-08-16	1	AUTOMOTIVE	0
1	3000889	2017-08-16	1	BABY CARE	0
2	3000890	2017-08-16	1	BEAUTY	2
Data type	int64	object	int64	object	int64
Unique values	28512	16	54	33	212
NaNs	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)

None

holidays:
(350, 6)

	date	type	locale	locale_name	description	transferred
0	2012-03-02	Holiday	Local	Manta	Fundacion de Manta	False
1	2012-04-01	Holiday	Regional	Cotopaxi	Provincializacion de Cotopaxi	False
2	2012-04-12	Holiday	Local	Cuenca	Fundacion de Cuenca	False
Data type	object	object	object	object	object	bool
Unique values	312	6	3	24	103	2
NaNs	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)

None

oil:
(1218, 2)

	date	dcoilwtico
0	2013-01-01	NaN
1	2013-01-02	93.14
2	2013-01-03	92.97
Data type	object	float64
Unique values	1218	998
NaNs	(0, 0.0 %)	(43, 3.5 %)

None

stores:
(54, 5)

	store_nbr	city	state	type	cluster
0	1	Quito	Pichincha	D	13
1	2	Quito	Pichincha	D	13
2	3	Quito	Pichincha	D	8
Data type	int64	object	object	object	int64
Unique values	54	22	16	5	17
NaNs	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)

None

transactions:
(83488, 3)

	date	store_nbr	transactions
0	2013-01-01	25	770
1	2013-01-02	1	2111
2	2013-01-02	2	2358
Data type	object	int64	int64
Unique values	1682	54	4993
NaNs	(0, 0.0 %)	(0, 0.0 %)	(0, 0.0 %)

None

Wykres dziedzin sprzedaży

families_mean = df.groupby('family').mean(numeric_only=True).sort_values('sales', ascending=False)
# families_mean

fig = px.bar(x=families_mean.index, y=families_mean['sales'])
fig.update_layout(template='plotly_dark', title='Średnia dzienna sprzedaż według rodzaju', yaxis_title='Średnia sprzedaż', xaxis_title='dziedzina')
check_png(fig)

Wykres ilości transakcji

fig = px.box(x=transactions['store_nbr'], y=transactions['transactions'])
fig.update_traces(jitter=0.75, marker=dict(size=3))
fig.update_layout(
    width=1500,
    template='plotly_dark', title='Wykres pudełkowy ilości transakcji w danym sklepie',
    xaxis_title="Numer sklepu",
    yaxis_title="Ilość transakcji")
check_png(fig)

Sklepy i dziedziny sprzedaży

# Kolory z pakietu plotly
all_colors = []

for i in px.colors.qualitative.Plotly:
    all_colors.append(i)
for i in px.colors.qualitative.Set1:
    all_colors.append(i)
for i in px.colors.qualitative.Set2:
    all_colors.append(i)
for i in px.colors.qualitative.Set3:
    all_colors.append(i)

stores = df.groupby(['store_nbr', 'family']).sum(numeric_only=True).reset_index()

color_map = dict(zip(stores['family'].unique(), all_colors[:33]))

# stores
# color_map

fig = go.Figure()
for family, color in color_map.items():
    subset_df = stores[stores['family'] == family]
    fig.add_trace(go.Bar(
        x=subset_df['store_nbr'],
        y=subset_df['sales'],
        text=subset_df['family'],
        marker_color=color,
        name=family
    ))

fig.update_layout(
    template='plotly_dark', width=1900, height=1200,
    title='Całkowita sprzedaż w danym sklepie, zależnie od dziedziny', xaxis_title="Numer sklepu", yaxis_title="Suma sprzedaży", legend_title='Dziedzina',  
    showlegend=True,  
    barmode='stack',
    xaxis=dict(
        tickmode='linear',  
        dtick=3  
    )
)
check_png(fig, PNG_PLOTS)

Promocje

Funckja tworząca trace dla każdej dziedziny sprzedaży i ilości promocji.

• trace to obiekt graficzny, zawierający dane do wizualizacji wykresu
• dane te zostają nałożone na "płótno"
• styl wykresu można określić przy tworzeniu trace-ów lub przy nakładaniu

Do wykreślenia sprzedaży i promocji wykorzystano:

• metodę subplots - pozwala na wyświetlenie wielu wykresów na jednej "figurze"
• parametr secondary_y - umożliwia nałożenie na jeden wykres dwie krzywe w różnych, osobnych skalach, w tym przypadku główna (lewa) oś Y to sprzedaż produktów danej dziedziny, a drugorzędna (prawa) oś Y to ilość produktów w promocji

prom_figs = []
def get_sales_and_prom(sales_data, family):
    q = sales_data.query(f'family == "{family}"')[['date', 'sales', 'onpromotion']].groupby('date').sum().reset_index()
    trace_p = go.Scatter(x=q['date'], y=q['onpromotion'], mode='markers', marker=dict(size=2))
    trace_s = go.Scatter(x=q['date'], y=q['sales'], line=dict(width=1))

    trace_set = [trace_s, trace_p, family]
    prom_figs.append(trace_set)

for i in df['family'].unique():
    get_sales_and_prom(df, i)

figs = prom_figs
fig = make_subplots(rows=4, cols=4, specs=[[{"secondary_y": True}]*4]*4)
counter = 0
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(figs[counter][0], row=i, col=j, )
        fig.add_trace(figs[counter][1], row=i, col=j, secondary_y=True)
        fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        counter += 1
print(counter)
fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False)
check_png(fig, PNG_PLOTS)

16

# niestety `secondary_y i szerokie `subplots` jakoś nie działało
figs = prom_figs
fig = make_subplots(rows=4, cols=4, specs=[[{"secondary_y": True}]*4]*4)
counter = 16
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(figs[counter][0], row=i, col=j, )
        fig.add_trace(figs[counter][1], row=i, col=j, secondary_y=True)
        fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        counter += 1

fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False)
check_png(fig, PNG_PLOTS)

fig = make_subplots(specs=[[{"secondary_y": True}]])
fig.add_trace(figs[32][0], )
fig.add_trace(figs[32][1], secondary_y=True)
fig.update_yaxes(title_text=figs[32][-1],)

fig.update_layout(height=400, width=800, template='plotly_dark', showlegend=False)
check_png(fig, PNG_PLOTS)

Wnioski z promocji:

• w większości brak wyraźnej korelacji

• najbardziej wyraźna w SCHOOL AND OFFICE SUPPLIES, w okolicach rozpoczęcia roku szkolnego (w Ekwadorze również początek września)

• LAWN AND GARDEN również wykazuje duży związek między wzrostem sprzedaży i promocjami

• podczas modelowania z lagami dla promocji wyniki były mniej trafne

Ropa / benzyna

• tutaj też zostaną wykorzystane drugorzędne osie Y
• prawdopodobnie takie rozpatrywanie nie jest konieczne, ale pozwoli na łatwe zidentyfikowanie, które dziedziny sprzedaży najbardziej korelują z ceną benzyny, tudzież ogólnym stanem gospodarki (gospodarka Ekwadoru silnie opiera się na wydobyciu ropy)

oil = oil.fillna(method='backfill')
oil_c = oil[:-12]
oil_c.tail()

	date	dcoilwtico
1201	2017-08-09	49.59
1202	2017-08-10	48.54
1203	2017-08-11	48.81
1204	2017-08-14	47.59
1205	2017-08-15	47.57

sales_per_day = df.groupby('date').sum(numeric_only=True)['sales']

fig = make_subplots(specs=[[{"secondary_y": True}]])
fig.add_trace(go.Scatter(x=sales_per_day.index, y=sales_per_day, name='Total sales per day'))
fig.add_trace(go.Scatter(x=oil['date'], y=oil['dcoilwtico'], name='Daily oil price'), secondary_y=True)
fig.update_yaxes(
    
)

fig.update_layout(height=800, width=1800, template='plotly_dark',)
check_png(fig, False)

figs = prom_figs
fig = make_subplots(rows=4, cols=4, specs=[[{"secondary_y": True}]*4]*4)
counter = 0
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(figs[counter][0], row=i, col=j, )
        fig.add_trace(go.Scatter(x=oil['date'], y=oil['dcoilwtico'],), row=i, col=j, secondary_y=True)
        fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        fig.update_yaxes(title_text='Oil', row=i, col=j, range=[-0, 120], secondary_y=True)
        counter += 1
print(counter)
fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False,
                  title='Sprzedaż i ropa/benzyna')
check_png(fig, PNG_PLOTS)

16

figs = prom_figs
fig = make_subplots(rows=4, cols=4, specs=[[{"secondary_y": True}]*4]*4)
counter = 16
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(figs[counter][0], row=i, col=j, )
        fig.add_trace(go.Scatter(x=oil['date'], y=oil['dcoilwtico'], name='Daily oil price'), row=i, col=j, secondary_y=True)
        fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        fig.update_yaxes(title_text='Oil', row=i, col=j, range=[-0, 120], secondary_y=True)
        counter += 1

fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False,
                  title='Sprzedaż i ropa/benzyna')

check_png(fig, PNG_PLOTS)

Korelacje

Korelacje dziedzin względem:

• promocji
• ropy/benzyny

corrs = df.groupby(['date','family']).sum(numeric_only=True)[['sales', 'onpromotion']].unstack().corr()
corrs = corrs.reset_index()
corrs['new_family'] = corrs['level_0'] + ' ' + corrs['family']
corrs = corrs.set_index('new_family')
corrs.drop(columns=['level_0', 'family'], level=0, inplace=True)

corrs = corrs.transpose().reset_index()
corrs['new_family'] = corrs['level_0'] + ' ' + corrs['family']
corrs
corrs = corrs.set_index('new_family')
corrs.drop(columns=['level_0', 'family'], inplace=True)

corrs.fillna(0, inplace=True)

# fig = px.imshow(corrs.round(decimals=2) * 100, # w skali 1-100 żeby zajmowało mniej miejsca w kwadracie
#                 text_auto=True,
#                 labels=dict(color="Correlation"),
#                 color_continuous_scale='RdBu_r', 
#                 color_continuous_midpoint=0)

# fig.update_traces(textfont_size=8)
# fig.update_xaxes(tickfont=dict(size=9))
# fig.update_yaxes(tickfont=dict(size=9))

# fig.update_layout(height=1200, width=1200,)

# check_png(fig)

grouped = df.groupby("family")

correlations = {}

for family, group in grouped:
    correlation = group["sales"].corr(group["onpromotion"])
    correlations[family] = correlation

fig = px.bar(pd.Series(correlations).sort_values(ascending=False))
fig.update_layout(template='plotly_dark', title='Korelacja sprzedaży względem promocji', 
                  width=1200, showlegend=False,
                  xaxis_title="Dziedzina sprzedaży",
                  yaxis_title="Korelacja",)
check_png(fig, PNG_PLOTS)

grouped = df.merge(right=oil_c).groupby("family")

correlations = {}

for family, group in grouped:
    correlation = group["sales"].corr(group["dcoilwtico"])
    correlations[family] = correlation

fig = px.bar(pd.Series(correlations).sort_values(ascending=True))
fig.update_layout(template='plotly_dark', title='Korelacja sprzedaży względem ceny beznyny', 
                  width=1200, showlegend=False,
                  xaxis_title="Dziedzina sprzedaży",
                  yaxis_title="Korelacja",)
check_png(fig, PNG_PLOTS)

Analiza czasowa

W analizie czasowej zostanie rozpatrzone jak sprzedaż zmienia się w czasie. Uwzględni dni tygodnia, średnią ruchomą, regresję liniową i wielomianową

def fix_dates(df, datename='date', set_date=True):
    df['date'] = pd.to_datetime(df[datename], format='%Y-%m-%d')
    df['date'] = df.date.dt.to_period('D')
    if set_date:
        df = df.set_index('date')   
    return df

grouped = df.copy()
grouped = fix_dates(grouped, 'date', False).groupby('date').mean(numeric_only=True)

grouped['day of week'] = grouped.index.dayofweek

Dni tygodnia według dokumentacji pandas

grouped

	id	store_nbr	sales	onpromotion	day of week
date
2013-01-01	890.5	27.5	1.409438	0.000000	1
2013-01-02	2672.5	27.5	278.390807	0.000000	2
2013-01-03	4454.5	27.5	202.840197	0.000000	3
2013-01-04	6236.5	27.5	198.911154	0.000000	4
2013-01-05	8018.5	27.5	267.873244	0.000000	5
...	...	...	...	...	...
2017-08-11	2992868.5	27.5	463.733851	7.956790	4
2017-08-12	2994650.5	27.5	444.798280	4.664422	5
2017-08-13	2996432.5	27.5	485.768618	5.209315	6
2017-08-14	2998214.5	27.5	427.004717	4.513468	0
2017-08-15	2999996.5	27.5	427.980884	5.951178	1

1684 rows × 5 columns

day_names = {0: 'Poniedziałek', 1: 'Wtorek', 2: 'Środa', 3: 'Czwartek', 4: 'Piątek', 5: 'Sobota', 6: 'Niedziela'}
grouped['day of week'] = grouped['day of week'].map(day_names)

grouped

	id	store_nbr	sales	onpromotion	day of week
date
2013-01-01	890.5	27.5	1.409438	0.000000	Wtorek
2013-01-02	2672.5	27.5	278.390807	0.000000	Środa
2013-01-03	4454.5	27.5	202.840197	0.000000	Czwartek
2013-01-04	6236.5	27.5	198.911154	0.000000	Piątek
2013-01-05	8018.5	27.5	267.873244	0.000000	Sobota
...	...	...	...	...	...
2017-08-11	2992868.5	27.5	463.733851	7.956790	Piątek
2017-08-12	2994650.5	27.5	444.798280	4.664422	Sobota
2017-08-13	2996432.5	27.5	485.768618	5.209315	Niedziela
2017-08-14	2998214.5	27.5	427.004717	4.513468	Poniedziałek
2017-08-15	2999996.5	27.5	427.980884	5.951178	Wtorek

1684 rows × 5 columns

fig = px.box(x=grouped['day of week'], y=grouped['sales'])
fig.update_traces(jitter=0.75, marker=dict(size=3))
fig.update_layout(
    width=1500,
    template='plotly_dark', title='Wykres pudełkowy sprzedaży zależnie od dnia tygodnia',
    xaxis_title="Dzień tygodnia",
    yaxis_title="Sprzedaż",
    xaxis=dict(
        categoryorder="array",
        categoryarray=list(day_names.values()))
    )
check_png(fig)

Wnioski ze sprzedaży zależnie od dnia tygodnia:

• największa sprzedaż w niedzielę i sobotę
• najmniejsza w czwartek
• zdecydowanie to rozróżnienie będzie istotne w modelowaniu

Trend

sales_per_day = sales_per_day.to_frame()
sales_per_day['time_step'] = np.arange(len(sales_per_day.index))
sales_per_day

	sales	time_step
date
2013-01-01	2511.618999	0
2013-01-02	496092.417944	1
2013-01-03	361461.231124	2
2013-01-04	354459.677093	3
2013-01-05	477350.121229	4
...	...	...
2017-08-11	826373.722022	1679
2017-08-12	792630.535079	1680
2017-08-13	865639.677471	1681
2017-08-14	760922.406081	1682
2017-08-15	762661.935939	1683

1684 rows × 2 columns

X = sales_per_day.loc[:, ['time_step']]
y = sales_per_day.loc[:, 'sales']

model = LinearRegression()
model.fit(X, y)

LinearRegression()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

y_pred_lr = pd.Series(model.predict(X), index=X.index)

all_holidays = holidays.query('locale == "National"').join(sales_per_day, on='date').fillna(0)

transfers = holidays.query('transferred == True')
transfers = transfers.join(sales_per_day, on='date').fillna(0)
transfers

	date	type	locale	locale_name	description	transferred	sales	time_step
19	2012-10-09	Holiday	National	Ecuador	Independencia de Guayaquil	True	0.000000	0.0
72	2013-10-09	Holiday	National	Ecuador	Independencia de Guayaquil	True	322529.418957	281.0
135	2014-10-09	Holiday	National	Ecuador	Independencia de Guayaquil	True	510519.923104	645.0
255	2016-05-24	Holiday	National	Ecuador	Batalla de Pichincha	True	606377.205216	1236.0
266	2016-07-25	Holiday	Local	Guayaquil	Fundacion de Guayaquil	True	697383.602092	1298.0
268	2016-08-10	Holiday	National	Ecuador	Primer Grito de Independencia	True	658457.436112	1314.0
297	2017-01-01	Holiday	National	Ecuador	Primer dia del ano	True	12082.500997	1457.0
303	2017-04-12	Holiday	Local	Cuenca	Fundacion de Cuenca	True	791762.312883	1558.0
312	2017-05-24	Holiday	National	Ecuador	Batalla de Pichincha	True	746303.627126	1600.0
324	2017-08-10	Holiday	National	Ecuador	Primer Grito de Independencia	True	651386.911970	1678.0
328	2017-09-28	Holiday	Local	Ibarra	Fundacion de Ibarra	True	0.000000	0.0
340	2017-12-06	Holiday	Local	Quito	Fundacion de Quito	True	0.000000	0.0

window_v = 30
SMA_30 = sales_per_day['sales'].rolling(window=window_v, 
                                     center=True, 
                                     min_periods=int(window_v/2)
                                     ).mean()
window_v = 365
SMA_365 = sales_per_day['sales'].rolling(window= window_v, 
                                     center=True, 
                                     min_periods=int(window_v/2)
                                     ).mean()

fig = go.Figure()
fig.add_trace(go.Scatter(x=sales_per_day.index, y=sales_per_day['sales'], mode='lines', name='Total sales', 
                         line=dict(width=1.5)))
fig.add_trace(go.Scatter(x=all_holidays['date'], y=all_holidays['sales'], mode='markers', name='All holidays'))
# fig.add_trace(go.Scatter(x=transfers['date'], y=transfers['sales'], mode='markers', name='Transferred holidays'))
fig.add_trace(go.Scatter(x=y_pred_lr.index, y=y_pred_lr, mode='lines', name='Regression fit', 
                         line=dict(width=3), marker=dict(color='#FECB52')))
fig.add_trace(go.Scatter(x=sales_per_day.index, y=SMA_30, 
                         name=f'Trend sprzedaży w przedziale czasowym 30 dni'))
fig.add_trace(go.Scatter(x=sales_per_day.index, y=SMA_365, 
                         name=f'Trend sprzedaży w przedziale czasowym 180 dni'))


fig.update_layout(template='plotly_dark', 
                  height=700, width=1800,
                  legend=dict(
                      orientation="h",
                      yanchor="bottom",
                      y=1.02,
                      xanchor="right",
                      x=1))

check_png(fig)

warnings.filterwarnings('ignore', category=UserWarning)

average_sales = sales_per_day['sales']
y = average_sales.copy()

dp = DeterministicProcess(average_sales.index, constant=False, order=3) # do 3 potęgi, krzywa wielomianowa zamiast linii 

X = dp.in_sample()

forecast_end = 32
X_fore = dp.out_of_sample(forecast_end) # kontynuacja krzywej na 16 kroków (dni)

average_sales.index

Index(['2013-01-01', '2013-01-02', '2013-01-03', '2013-01-04', '2013-01-05',
       '2013-01-06', '2013-01-07', '2013-01-08', '2013-01-09', '2013-01-10',
       ...
       '2017-08-06', '2017-08-07', '2017-08-08', '2017-08-09', '2017-08-10',
       '2017-08-11', '2017-08-12', '2017-08-13', '2017-08-14', '2017-08-15'],
      dtype='object', name='date', length=1684)

X_fore_steps = X_fore.index.to_list()
print(X_fore_steps)

[1685, 1686, 1687, 1688, 1689, 1690, 1691, 1692, 1693, 1694, 1695, 1696, 1697, 1698, 1699, 1700, 1701, 1702, 1703, 1704, 1705, 1706, 1707, 1708, 1709, 1710, 1711, 1712, 1713, 1714, 1715, 1716]

start_date = datetime.strptime("2017-08-16", "%Y-%m-%d")
end_date = start_date + timedelta(days = forecast_end) 

datelist = pd.date_range(start="2017-08-16", end=end_date)
datelist

DatetimeIndex(['2017-08-16', '2017-08-17', '2017-08-18', '2017-08-19',
               '2017-08-20', '2017-08-21', '2017-08-22', '2017-08-23',
               '2017-08-24', '2017-08-25', '2017-08-26', '2017-08-27',
               '2017-08-28', '2017-08-29', '2017-08-30', '2017-08-31',
               '2017-09-01', '2017-09-02', '2017-09-03', '2017-09-04',
               '2017-09-05', '2017-09-06', '2017-09-07', '2017-09-08',
               '2017-09-09', '2017-09-10', '2017-09-11', '2017-09-12',
               '2017-09-13', '2017-09-14', '2017-09-15', '2017-09-16',
               '2017-09-17'],
              dtype='datetime64[ns]', freq='D')

X_fore['date'] = datelist[:forecast_end]
X_fore = X_fore.set_index('date')
X_fore

	trend	trend_squared	trend_cubed
date
2017-08-16	1685.0	2839225.0	4.784094e+09
2017-08-17	1686.0	2842596.0	4.792617e+09
2017-08-18	1687.0	2845969.0	4.801150e+09
2017-08-19	1688.0	2849344.0	4.809693e+09
2017-08-20	1689.0	2852721.0	4.818246e+09
2017-08-21	1690.0	2856100.0	4.826809e+09
2017-08-22	1691.0	2859481.0	4.835382e+09
2017-08-23	1692.0	2862864.0	4.843966e+09
2017-08-24	1693.0	2866249.0	4.852560e+09
2017-08-25	1694.0	2869636.0	4.861163e+09
2017-08-26	1695.0	2873025.0	4.869777e+09
2017-08-27	1696.0	2876416.0	4.878402e+09
2017-08-28	1697.0	2879809.0	4.887036e+09
2017-08-29	1698.0	2883204.0	4.895680e+09
2017-08-30	1699.0	2886601.0	4.904335e+09
2017-08-31	1700.0	2890000.0	4.913000e+09
2017-09-01	1701.0	2893401.0	4.921675e+09
2017-09-02	1702.0	2896804.0	4.930360e+09
2017-09-03	1703.0	2900209.0	4.939056e+09
2017-09-04	1704.0	2903616.0	4.947762e+09
2017-09-05	1705.0	2907025.0	4.956478e+09
2017-09-06	1706.0	2910436.0	4.965204e+09
2017-09-07	1707.0	2913849.0	4.973940e+09
2017-09-08	1708.0	2917264.0	4.982687e+09
2017-09-09	1709.0	2920681.0	4.991444e+09
2017-09-10	1710.0	2924100.0	5.000211e+09
2017-09-11	1711.0	2927521.0	5.008988e+09
2017-09-12	1712.0	2930944.0	5.017776e+09
2017-09-13	1713.0	2934369.0	5.026574e+09
2017-09-14	1714.0	2937796.0	5.035382e+09
2017-09-15	1715.0	2941225.0	5.044201e+09
2017-09-16	1716.0	2944656.0	5.053030e+09

model = LinearRegression()
model.fit(X, y)

y_pred = pd.Series(model.predict(X), index=X.index)
y_fore = pd.Series(model.predict(X_fore), index=X_fore.index)

window_v = 30
SMA_30 = sales_per_day['sales'].rolling(window=window_v, 
                                     center=True, 
                                     min_periods=int(window_v/2)
                                     ).mean()
window_v = 180
SMA_365 = sales_per_day['sales'].rolling(window= window_v, 
                                     center=True, 
                                     min_periods=int(window_v/2)
                                     ).mean()

fig = go.Figure()
fig.add_trace(go.Scatter(x=sales_per_day.index, y=sales_per_day['sales'], mode='lines', name='Total sales', 
                         line=dict(width=1.5)))
fig.add_trace(go.Scatter(x=all_holidays['date'], y=all_holidays['sales'], mode='markers', name='All holidays'))
fig.add_trace(go.Scatter(x=transfers['date'], y=transfers['sales'], mode='markers', name='Transferred holidays'))
fig.add_trace(go.Scatter(x=y_pred_lr.index, y=y_pred_lr, mode='lines', name='Regression fit', 
                         line=dict(width=2), marker=dict(color='#b4913e')))
fig.add_trace(go.Scatter(x=y_pred.index, y=y_pred, mode='lines', name='Polynomial-regression fit', 
                         line=dict(width=2), marker=dict(color='#FECB52')))
fig.add_trace(go.Scatter(x=sales_per_day.index, y=SMA_30, 
                         name=f'Trend sprzedaży w przedziale czasowym 30 dni'))
fig.add_trace(go.Scatter(x=sales_per_day.index, y=SMA_365, 
                         name=f'Trend sprzedaży w przedziale czasowym 180 dni'))
fig.add_trace(go.Scatter(x=y_fore.index, y=y_fore, 
                         name=f'Prognoza sprzedaży na {forecast_end} dni'))

fig.update_xaxes(range=[datetime.strptime("2012-12-01", "%Y-%m-%d"),
                        end_date])

fig.update_layout(template='plotly_dark', 
                  height=700, width=1800,
                  legend=dict(
                      orientation="h",
                      yanchor="bottom",
                      y=1.02,
                      xanchor="right",
                      x=1))

check_png(fig, PNG_PLOTS)

Poniżej znajdują się wykreślone zestawy innych dziedzin sprzedaży

q_test = df.groupby(['date', 'family']).sum().query('family == "AUTOMOTIVE"')['sales'].to_frame()
q_test['time_step'] = np.arange(len(q_test.index))
q_test

		sales	time_step
date	family
2013-01-01	AUTOMOTIVE	0.0	0
2013-01-02	AUTOMOTIVE	255.0	1
2013-01-03	AUTOMOTIVE	161.0	2
2013-01-04	AUTOMOTIVE	169.0	3
2013-01-05	AUTOMOTIVE	342.0	4
...	...	...	...
2017-08-11	AUTOMOTIVE	441.0	1679
2017-08-12	AUTOMOTIVE	403.0	1680
2017-08-13	AUTOMOTIVE	481.0	1681
2017-08-14	AUTOMOTIVE	292.0	1682
2017-08-15	AUTOMOTIVE	337.0	1683

1684 rows × 2 columns

def get_lr(query_part):
    q_test = query_part.to_frame()
    q_test['time_step'] = np.arange(len(q_test.index))

    X = q_test.loc[:, ['time_step']]  # features
    y = q_test.loc[:, 'sales']  # target

    model = LinearRegression()
    model.fit(X, y)

    y_pred_lr = pd.Series(model.predict(X), index=X.index)
    # y_pred_lr = y_pred_lr.to_frame().reset_index().drop(columns='family').set_index('date')[0]

    return y_pred_lr

get_lr(df.groupby(['date', 'family']).sum().query('family == "AUTOMOTIVE"')['sales'])

date        family    
2013-01-01  AUTOMOTIVE    241.322651
2013-01-02  AUTOMOTIVE    241.427398
2013-01-03  AUTOMOTIVE    241.532144
2013-01-04  AUTOMOTIVE    241.636891
2013-01-05  AUTOMOTIVE    241.741637
                             ...    
2017-08-11  AUTOMOTIVE    417.191855
2017-08-12  AUTOMOTIVE    417.296601
2017-08-13  AUTOMOTIVE    417.401348
2017-08-14  AUTOMOTIVE    417.506094
2017-08-15  AUTOMOTIVE    417.610840
Length: 1684, dtype: float64

def get_SMA(window, sales):
    SMA = sales.rolling(window=window, 
                                     center=True, 
                                     min_periods=int(window/2)
                                     ).mean()
    return SMA

get_SMA(30, df.query(f"family == 'AUTOMOTIVE'").groupby('date').sum(numeric_only=True)['sales'])

date
2013-01-01    208.533333
2013-01-02    206.187500
2013-01-03    205.352941
2013-01-04    203.388889
2013-01-05    207.473684
                 ...    
2017-08-11    406.700000
2017-08-12    410.263158
2017-08-13    413.944444
2017-08-14    406.470588
2017-08-15    398.812500
Name: sales, Length: 1684, dtype: float64

def get_dt_forecast(start_date=start_date, end_date=end_date, datelist=datelist, sales=None):
    y = sales.copy()  # the target
    dp = DeterministicProcess(average_sales.index, constant=False, order=3)
    X = dp.in_sample()
    X_fore = dp.out_of_sample(forecast_end)
    X_fore_steps = X_fore.index.to_list()

    X_fore['date'] = datelist[:forecast_end]
    X_fore = X_fore.set_index('date')

    model = LinearRegression()
    model.fit(X, y)

    y_pred = pd.Series(model.predict(X), index=X.index)
    y_fore = pd.Series(model.predict(X_fore), index=X_fore.index)

    return [y_pred, y_fore]

# get_dt_forecast(30, sales = df.query(f"family == 'AUTOMOTIVE'").groupby('date').sum(numeric_only=True)['sales']) 

figs = []
for i in df['family'].unique():
    q = df.query(f"family == '{i}'").groupby('date').sum(numeric_only=True)['sales']
    trace_base = go.Scatter(x=q.index, y=q, name=f'Sprzedaż {i}', line=dict(width=1), marker=dict(color=px.colors.qualitative.Plotly[0]))

    lr = get_lr(q)
    trace_lr = go.Scatter(x=lr.index, y=lr, name=f'Regresja liniowa {i}', line=dict(width=2), marker=dict(color='#b4913e'))

    dt_forecast = get_dt_forecast(sales=q)
    trace_dt = go.Scatter(x=dt_forecast[0].index, y=dt_forecast[0], name=f'Trend {i}', line=dict(width=1.5), marker=dict(color='#FECB52'))
    trace_fr = go.Scatter(x=dt_forecast[1].index, y=dt_forecast[1], name=f'Prognoza {i}', line=dict(width=1.5), marker=dict(color=px.colors.qualitative.Plotly[2]))

    c_SMA_30 = get_SMA(30, q)
    c_SMA_180 = get_SMA(180, q)
    trace_SMA_30 = go.Scatter(x=c_SMA_30.index, y=c_SMA_30, name=f'Trend sprzedaży w przedziale czasowym 30 dni', marker=dict(color=px.colors.qualitative.Plotly[6]))
    trace_SMA_180 = go.Scatter(x=c_SMA_180.index, y=c_SMA_180, name=f'Trend sprzedaży w przedziale czasowym 180 dni', marker=dict(color=px.colors.qualitative.Plotly[4]))



    trace_set = []
    trace_set.append(trace_base)
    trace_set.append(trace_lr)
    trace_set.append(trace_dt)
    trace_set.append(trace_fr)
    trace_set.append(trace_SMA_30)
    trace_set.append(trace_SMA_180)

    trace_set.append(i)
    figs.append(trace_set)

fig = make_subplots(rows=4, cols=4)
counter = 0
for i in range(1,5):
    for j in range(1,5):
        for k in range(0,6):
            fig.add_trace(figs[counter][k], row=i, col=j, )
            fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        counter += 1


fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False)

                  

check_png(fig, PNG_PLOTS)

fig = make_subplots(
  rows=5, cols=4,
  specs=[
      [{}, {}, {}, {}],
      [{}, {}, {}, {}],
      [{}, {}, {}, {}],
      [{}, {}, {"rowspan":1, "colspan": 2}, None,],
      [{"rowspan":1, "colspan": 2}, None, {"rowspan":1, "colspan": 2}, None]
         ],
  print_grid=False)

counter = 16
for i in range(1,4):
    for j in range(1,5):
        for k in range(0,6):
            fig.add_trace(figs[counter][k], row=i, col=j, )
            fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        counter += 1

for j in range(1,4):
    for k in range(0,6):
        fig.add_trace(figs[counter][k], row=4, col=j, )
        fig.update_yaxes(title_text=figs[counter][-1], row=4, col=j, )
    counter += 1

for k in range(0,6):
        fig.add_trace(figs[counter][k], row=5, col=1, )
        fig.update_yaxes(title_text=figs[counter][-1], row=5, col=1, )
counter += 1

for k in range(0,6):
        fig.add_trace(figs[counter][k], row=5, col=3, )
        fig.update_yaxes(title_text=figs[counter][-1], row=5, col=3, )

fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False)
check_png(fig, PNG_PLOTS)

Wnioski z trendu:

• większość sprzedaży stopniowo rośnie z czasem
• koniec 2017 roku często wskazuje na spadek
• w wielu przypadkach mamy zupełne braki sprzedaży, które wprowadzają błędy do analizy - modelowanie zostanie przeprowadzona na podzbiorach przefiltrowanych do opdowiednich dat

Dziedziny bez zarzutu:

['AUTOMOTIVE', 'BEAUTY', 'BREAD/BAKERY', 'CLEANING', 'DELI', 'EGGS', 'GROCERY I', 'GROCERY II', 'HARDWARE',]

['HOME APPLIANCES', 'LINGERIE', 'PERSONAL CARE', 'PREPARED FOODS', 'SEAFOOD']

Dziedziny dobre, ale do ucięcia w czasie:

• BABY CARE - od 2014-03-1 i edytować wartość z 2015-08-4 lub lepiej zacząć od 2015-12-7, sprawdzić onpromotion
• BEVERAGES - od 2015-05-28
• CELEBRATION - od 2015-06-01
• DIARY - od 2013-11-7
• FROZEN FOODS - może normalnie, sprawdzić dokąd sięgają dane zestawu testowego, może zrobić one-hot encoding dla grudnia?
• HOME AND KITCHEN I - od 2014-09-1
• HOME AND KITCHEN II - od 2015-10-31
• HOME CARE - od 2015-05-4
• LADIESWEAR - od 2015-06-1
• LAWN AND GARDEN - od 2016-12-3, może zmienić wartości dla 2017-02-13 -14 i 2017-05-12 -14
• LIQUOR,WINE,BEER - od 2016-05-10 - wcześniej alkohol nie mógł byś sprzedawany w niedziele, od mniej więcej tego momentu już może być sprzedawany w każdy dzień tygodnia
• MAGAZINES - od 2015-09-30
• MEATS - edytować dla 2016-10-7 lub sprawdzić onpromotion
• PET SUPPLIES - od 2015-06-1
• PLAYERS AND ELECTRONICS - od 2015-06-1
• POULTRY - od 2013-11-4
• PRODUCE - od 2015-06-1

Słabe / niepewne dziedziny:

• BOOKS - może od 2017-01-1
• SCHOOLS AND OFFICE SUPPLIES - może wprowadzić wydarzenie początek roku szkolnego i może lag dla niego?

Trend - ropa/benzyna

oil.head()

	date	dcoilwtico
0	2013-01-01	93.14
1	2013-01-02	93.14
2	2013-01-03	92.97
3	2013-01-04	93.12
4	2013-01-07	93.20

oil_SMA30 = get_SMA(30, oil.set_index('date')['dcoilwtico'])
oil_SMA30_trace = go.Scatter(x=oil_SMA30.index, y=oil_SMA30, name ='Trend ceny ropy')

fig = make_subplots(rows=4, cols=4, specs=[[{"secondary_y": True}]*4]*4)
counter = 0
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(figs[counter][4], row=i, col=j, )
        fig.add_trace(oil_SMA30_trace, row=i, col=j, secondary_y=True)
        fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        fig.update_yaxes(title_text="OIL SMA", row=i, col=j, secondary_y=True)
        counter += 1

fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False,
                  title='Średnia ruchoma dla 30 dni')

check_png(fig, PNG_PLOTS)

fig = make_subplots(rows=4, cols=4, specs=[[{"secondary_y": True}]*4]*4)
counter = 16
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(figs[counter][4], row=i, col=j, )
        fig.add_trace(oil_SMA30_trace, row=i, col=j, secondary_y=True)
        fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        fig.update_yaxes(title_text="OIL SMA", row=i, col=j, secondary_y=True)
        counter += 1

fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False,
                  title='Średnia ruchoma dla 30 dni',)

check_png(fig, PNG_PLOTS)

fig = make_subplots(specs=[[{"secondary_y": True}]])
fig.add_trace(figs[32][4],)
fig.add_trace(oil_SMA30_trace, secondary_y=True, )
fig.update_yaxes(title_text=figs[32][-1],)
fig.update_yaxes(title_text="OIL SMA", secondary_y=True)

fig.update_layout(template='plotly_dark', 
                  height=400, width=1200,
                  legend=dict(
                      orientation="h",
                      yanchor="bottom",
                      y=1.02,
                      xanchor="right",
                      x=1))
check_png(fig, PNG_PLOTS)

Sezonowość

• dzień tygodnia zostanie zakodowany one-hot (indicators)
• sezonowość zostanie zbadana poprzez periodogram

def trace_periodogram(ts, detrend='linear', name='Periodogram'):
    from scipy.signal import periodogram
    fs = pd.Timedelta("365D") / pd.Timedelta("1D")
    freqencies, spectrum = periodogram(
        ts,
        fs=fs,
        detrend=detrend,
        window="boxcar",
        scaling='spectrum',
    )
    df = pd.DataFrame({'Frequency': freqencies, 'Spectrum': spectrum})
    
    trace = go.Scatter(
        x=df['Frequency'],
        y=df['Spectrum'],
        mode='lines',
        name=name,
    )
    return trace

periodograms = []
for i in df['family'].unique():
    q = df.query(f'family == "{i}"').groupby('date').sum(numeric_only=True)['sales']
    periodograms.append([trace_periodogram(q, name=i), i])

periodograms = [[trace, 'Var: ' + title] for trace, title in periodograms]

fig = make_subplots(rows=4, cols=4)
counter_a = 0
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(periodograms[counter_a][0], row=i, col=j, )
        fig.update_yaxes(title_text=periodograms[counter_a][-1], row=i, col=j, )  
        counter_a += 1

fig.update_xaxes(tickangle=30, tickvals=[1, 2, 4, 6, 12, 26, 52, 104],
        ticktext=[
            "Annual (1)",
            "Semiannual (2)",
            "Quarterly (4)",
            "Bimonthly (6)",
            "Monthly (12)",
            "Biweekly (26)",
            "Weekly (52)",
            "Semiweekly (104)",],)

fig.update_xaxes(type="log")
fig.update_layout(width=1500)
fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False,
                  title='Sezonowość - wariancja sprzedaży w zależności od przedziału czasowego')
check_png(fig, PNG_PLOTS)

fig = make_subplots(
  rows=5, cols=4,
  specs=[
      [{}, {}, {}, {}],
      [{}, {}, {}, {}],
      [{}, {}, {}, {}],
      [{}, {}, {"rowspan":1, "colspan": 2}, None,],
      [{"rowspan":1, "colspan": 2}, None, {"rowspan":1, "colspan": 2}, None]
         ],

  print_grid=False)

counter_b = 16

for i in range(1,4):
    for j in range(1,5):
        fig.add_trace(periodograms[counter_b][0], row=i, col=j, )
        fig.update_yaxes(title_text=periodograms[counter_b][-1], row=i, col=j, )  
        counter_b += 1

for j in range(1,4):
    fig.add_trace(periodograms[counter_b][0], row=4, col=j, )
    fig.update_yaxes(title_text=periodograms[counter_b][-1], row=4, col=j, )  
    counter_b += 1

fig.add_trace(periodograms[counter_b][0], row=5, col=1, )
fig.update_yaxes(title_text=periodograms[counter_b][-1], row=5, col=1, )  
counter_b += 1
fig.add_trace(periodograms[counter_b][0], row=5, col=3, )
fig.update_yaxes(title_text=periodograms[counter_b][-1], row=5, col=3, )  

fig.update_xaxes(
        tickvals=[1, 2, 4, 6, 12, 26, 52, 104],
        ticktext=[
            "Rocznie (1)",
            "Półrocznie (2)",
            "Kwartalnie (4)",
            "Co 2 miesiące (6)",
            "Miesięcznie (12)",
            "Co 2 tygodnie (26)",
            "Tygodniowo (52)",
            "x2 w tygodniu (104)",
        ],
        tickangle=30,
    )
fig.update_xaxes(type="log")
fig.update_layout(width=1500)

fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False,
                  title='Sezonowość - wariancja sprzedaży w zależności od przedziału czasowego')
check_png(fig, PNG_PLOTS)

Sezonowość wnioski

(wnioski również oparte na ogólnej historii sprzedaży)

• wiekszość najczęściej wybieranych produktów, czyli najczęściej artykuły spożywcze - GROCERY I, MEATS, SEAFOOD etc. mają bardzo dobrze widoczną sezonowość tygodniową

  • duża wariancja wewnątrz tygodnia - najmniejsze zakupy w czwartki, największe w niedziele jak już wcześniej stwierdzono
  • tu będzie one-hot encoding dni tygodnia czyli tzw. time indicators

• dziedziny w których widać największą sezonowość w dużych okresach czasu to:

  • FROZEN FOOD - ogromny skok przed Świętami Bożego Narodzenia
  • SCHOOL AND OFFICE SUPPLIES - tu głównie wrzesień - rok szkolny zaczyna się we wrześniu, ale mamy też wzrost sprzedaży w maju/kwietniu 2016 i 2017 roku

• duża sezonowość w skali ponad roku wynika też z okresów gdzie sprzedaż nie odbywała się, lub sprzedaż tych produktów nie była rejestrowana albo była inaczej kwalifikowana

  • przykładowo BABY CARE ma niemal zerową sprzedaż do początku 2014 roku, oraz zupełny brak sprzedaż od stycznia do maja 2015 roku
  • podobnie PET SUPPLIES

• przykładem innego przypadku może być drób: POULTRY - do końca października 2014 roku ilość sprzedaży była około dwa razy mniejsza niż po październiku 2014 roku, przypuszczam że taki gwałtowny wzrost mógłby być spowodowany przez:

  • rozpoczęcie sprzedaży w nowych sklepach / zmianę typu sklepu
  • zmianę w cenie produktów, które np. dotychczas miały zaniżoną cenę lub zaszła zmiana w metodzie ich opodatkowania - drób i produkty mleczne są "codziennymi" typami produktów, bez których większość ludzi nie może się obejść (w szczególności drób jest w większości tańszym mięsem)
  • zmianą strategii firmy - np. dana firma, z której pochodzi ten zestaw danych, zaczęła sprowadzać więcej tych produktów do Ekwadoru, przez co ich sprzedaż wzrosła

Przygotowanie zestawów treningowych

• ta cześć notesu wykonywana w osobno od eksploracji danych
• idea przygotowania polega na rozdzieleniu pierwotnego zestawu testowego na 54 podzbiory dla każdego ze sklepów
• przed podziałem utworzono dodatkowe cechy, takie jak dzień tygodnia, miesiąc, rok, dzień miesiąca itd. oraz trend i cykle Fouriera
• oil oraz święta będą dołączane przy każdym podziale podzbioru na dziedziny i prognozowaniu

Innymi słowy, prognozowanie będzie odbywało się dla każdej osobnej dziedziny w każdym osobnym sklepie, czyli zostanie wykonane $33 \times 54 = 1782$ osobnych prognoz

# czy przeprowadzić zapis w tym pokazowym zeszycie
zapis = False

Holiday dummies

• dummies świąt to zero-jedynkowo zapisana tabela, w której każde święto stanowi osobną kolumnę, a wiersze reprezentowane są przez daty
• święta i wydarzenia zostały odfiltrowane do tych, które nie zostały przeniesione i tych które są obchodzone w całym kraju

train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')
oil = pd.read_csv('oil.csv')
holidays = pd.read_csv('holidays_events.csv')

train = fix_dates(train)
test = fix_dates(test)
oil = fix_dates(oil)
holidays = fix_dates(holidays)

oil

	dcoilwtico
date
2013-01-01	NaN
2013-01-02	93.14
2013-01-03	92.97
2013-01-04	93.12
2013-01-07	93.20
...	...
2017-08-25	47.65
2017-08-28	46.40
2017-08-29	46.46
2017-08-30	45.96
2017-08-31	47.26

1218 rows × 1 columns

holidays

	type	locale	locale_name	description	transferred
date
2012-03-02	Holiday	Local	Manta	Fundacion de Manta	False
2012-04-01	Holiday	Regional	Cotopaxi	Provincializacion de Cotopaxi	False
2012-04-12	Holiday	Local	Cuenca	Fundacion de Cuenca	False
2012-04-14	Holiday	Local	Libertad	Cantonizacion de Libertad	False
2012-04-21	Holiday	Local	Riobamba	Cantonizacion de Riobamba	False
...	...	...	...	...	...
2017-12-22	Additional	National	Ecuador	Navidad-3	False
2017-12-23	Additional	National	Ecuador	Navidad-2	False
2017-12-24	Additional	National	Ecuador	Navidad-1	False
2017-12-25	Holiday	National	Ecuador	Navidad	False
2017-12-26	Additional	National	Ecuador	Navidad+1	False

350 rows × 5 columns

dum_holidays = pd.get_dummies(holidays.query('transferred == @False & locale == "National"')[['description']], drop_first=False).astype('uint8')

dum_holidays.head(3)

	description_Batalla de Pichincha	description_Black Friday	description_Carnaval	description_Cyber Monday	description_Dia de Difuntos	description_Dia de la Madre	description_Dia de la Madre-1	description_Dia del Trabajo	description_Inauguracion Mundial de futbol Brasil	description_Independencia de Cuenca	...	description_Terremoto Manabi+5	description_Terremoto Manabi+6	description_Terremoto Manabi+7	description_Terremoto Manabi+8	description_Terremoto Manabi+9	description_Traslado Batalla de Pichincha	description_Traslado Independencia de Guayaquil	description_Traslado Primer Grito de Independencia	description_Traslado Primer dia del ano	description_Viernes Santo
date
2012-08-10	0	0	0	0	0	0	0	0	0	0	...	0	0	0	0	0	0	0	0	0	0
2012-10-12	0	0	0	0	0	0	0	0	0	0	...	0	0	0	0	0	0	1	0	0	0
2012-11-02	0	0	0	0	1	0	0	0	0	0	...	0	0	0	0	0	0	0	0	0	0

3 rows × 72 columns

if zapis:
    dum_holidays.to_csv('holiday_dummies.csv', index=True)

Feature engineering

• tu zostaną dodane wcześniej wymienione cechy

Dzień tygodnia

train['day_of_week'] = train.index.dayofweek
day_names = {0: 'Monday', 1: 'Tuesday', 2: 'Wednesday', 3: 'Thursday', 4: 'Friday', 5: 'Saturday', 6: 'Sunday'}
train['day_of_week'] = train['day_of_week'].map(day_names)

dum_day = pd.get_dummies(train["day_of_week"], drop_first=False)
train = pd.concat([train, dum_day], axis=1)
train.drop("day_of_week", axis=1, inplace=True)

test['day_of_week'] = test.index.dayofweek
day_names = {0: 'Monday', 1: 'Tuesday', 2: 'Wednesday', 3: 'Thursday', 4: 'Friday', 5: 'Saturday', 6: 'Sunday'}
test['day_of_week'] = test['day_of_week'].map(day_names)

dum_day = pd.get_dummies(test["day_of_week"], drop_first=False)
test = pd.concat([test, dum_day], axis=1)
test.drop("day_of_week", axis=1, inplace=True)

Inne okresy

• kwartał, miesiąc, rok, dzień w roku, etc.

def get_date_details(df):
    df['quarter'] = df.index.quarter
    df['month'] = df.index.month
    df['year'] = df.index.year
    df['dayofyear'] = df.index.dayofyear
    df['dayofmonth'] = df.index.day
    df['weekofyear'] = df.index.week
    return df

train = get_date_details(train)
test = get_date_details(test)

test.head()

	id	store_nbr	family	onpromotion	Friday	Monday	Saturday	Sunday	Thursday	Tuesday	Wednesday	quarter	month	year	dayofyear	dayofmonth	weekofyear
date
2017-08-16	3000888	1	AUTOMOTIVE	0	0	0	0	0	0	0	1	3	8	2017	228	16	33
2017-08-16	3000889	1	BABY CARE	0	0	0	0	0	0	0	1	3	8	2017	228	16	33
2017-08-16	3000890	1	BEAUTY	2	0	0	0	0	0	0	1	3	8	2017	228	16	33
2017-08-16	3000891	1	BEVERAGES	20	0	0	0	0	0	0	1	3	8	2017	228	16	33
2017-08-16	3000892	1	BOOKS	0	0	0	0	0	0	0	1	3	8	2017	228	16	33

Ropa/benzyna

oil.fillna(method='bfill', inplace=True) # w niektóre dni brakuje informacji o cenie - są NaN

train = train.join(oil.fillna(method='bfill'), how='left').fillna(method='bfill') # dołączenie tabelki do zestawów, uzupełnienie brakujących danych - nie ma niektórych dni w tabeli `oil`
test = test.join(oil.fillna(method='bfill'), how='left').fillna(method='bfill')

Cykle Fouriera

• cykle fouriera pozwalają na wielokrotnie oszczędniejszy zapis cykliczności niż one-hot encoding każdej daty / kroku
• opiera się na parach krzywych sinusoid i cosinusoid o określonej "długości", odwzorowywujące sezonowość danych wartości
• ta część również była wykonywana w osobnym notesie

Optymalne typy danych przed zapisem podzbiorów

train['id'] = train['id'].astype('uint32')
train['store_nbr'] = train['store_nbr'].astype('uint8')
train['family'] = train['family'].astype('category')
train['onpromotion'] = train['onpromotion'].astype('uint16')

train['quarter'] = train['quarter'].astype('uint8')
train['month'] = train['month'].astype('uint8')
train['year'] = train['year'].astype('uint16')
train['dayofyear'] = train['dayofyear'].astype('uint16')
train['dayofmonth'] = train['dayofmonth'].astype('uint8')
train['weekofyear'] = train['weekofyear'].astype('uint8')

test['id'] = test['id'].astype('uint32')
test['store_nbr'] = test['store_nbr'].astype('uint8')
test['family'] = test['family'].astype('category')
test['onpromotion'] = test['onpromotion'].astype('uint16')

test['quarter'] = test['quarter'].astype('uint8')
test['month'] = test['month'].astype('uint8')
test['year'] = test['year'].astype('uint16')
test['dayofyear'] = test['dayofyear'].astype('uint16')
test['dayofmonth'] = test['dayofmonth'].astype('uint8')
test['weekofyear'] = test['weekofyear'].astype('uint8')

fourier = CalendarFourier(freq="A", order=10)  # 10 par sin/cos dla corocznej sezonowości ("A" - Annual)
dp = DeterministicProcess(
    index=train.index,
    constant=True,               # cecha dla wyrazu wolnego
    order=1,                     # trend (order=1 oznacza regresję liniową)
    seasonal=True,               # sezonowość tygodniowa (indicators)
    additional_terms=[fourier],  # sezonowość w ciągu roku (fourier)
    drop=True,                   
)

fourier_in_sample = dp.in_sample()  
fourier_in_sample

	const	trend	s(2,7)	s(3,7)	s(4,7)	s(5,7)	s(6,7)	s(7,7)	sin(1,freq=A-DEC)	cos(1,freq=A-DEC)	...	sin(6,freq=A-DEC)	cos(6,freq=A-DEC)	sin(7,freq=A-DEC)	cos(7,freq=A-DEC)	sin(8,freq=A-DEC)	cos(8,freq=A-DEC)	sin(9,freq=A-DEC)	cos(9,freq=A-DEC)	sin(10,freq=A-DEC)	cos(10,freq=A-DEC)
date
2013-01-01	1.0	1.0	0.0	0.0	0.0	0.0	0.0	0.0	0.000000	1.000000	...	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000
2013-01-01	1.0	2.0	1.0	0.0	0.0	0.0	0.0	0.0	0.000000	1.000000	...	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000
2013-01-01	1.0	3.0	0.0	1.0	0.0	0.0	0.0	0.0	0.000000	1.000000	...	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000
2013-01-01	1.0	4.0	0.0	0.0	1.0	0.0	0.0	0.0	0.000000	1.000000	...	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000
2013-01-01	1.0	5.0	0.0	0.0	0.0	1.0	0.0	0.0	0.000000	1.000000	...	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
2017-08-15	1.0	3000884.0	0.0	0.0	0.0	1.0	0.0	0.0	-0.680773	-0.732494	...	-0.976011	-0.217723	0.863142	-0.504961	-0.288482	0.957485	-0.440519	-0.897743	0.933837	0.357698
2017-08-15	1.0	3000885.0	0.0	0.0	0.0	0.0	1.0	0.0	-0.680773	-0.732494	...	-0.976011	-0.217723	0.863142	-0.504961	-0.288482	0.957485	-0.440519	-0.897743	0.933837	0.357698
2017-08-15	1.0	3000886.0	0.0	0.0	0.0	0.0	0.0	1.0	-0.680773	-0.732494	...	-0.976011	-0.217723	0.863142	-0.504961	-0.288482	0.957485	-0.440519	-0.897743	0.933837	0.357698
2017-08-15	1.0	3000887.0	0.0	0.0	0.0	0.0	0.0	0.0	-0.680773	-0.732494	...	-0.976011	-0.217723	0.863142	-0.504961	-0.288482	0.957485	-0.440519	-0.897743	0.933837	0.357698
2017-08-15	1.0	3000888.0	1.0	0.0	0.0	0.0	0.0	0.0	-0.680773	-0.732494	...	-0.976011	-0.217723	0.863142	-0.504961	-0.288482	0.957485	-0.440519	-0.897743	0.933837	0.357698

3000888 rows × 28 columns

Podział na podzbiory

train = pd.read_csv('train w features.csv',
                 dtype = {'id': 'uint32', 
                          'store_nbr': 'uint8', 
                          'family': 'category',
                          'onpromotion': 'uint16',

                          'quarter': 'uint8', 
                          'month': 'uint8', 
                          'year': 'uint16',
                          'dayofyear': 'uint16',
                          'dayofmonth': 'uint8',
                          'weekofyear': 'uint8',

                          'Friday' : 'uint8', 
                          'Monday' : 'uint8', 
                          'Saturday' : 'uint8', 
                          'Sunday' : 'uint8', 
                          'Thursday' : 'uint8', 
                          'Tuesday' : 'uint8', 
                          'Wednesday' : 'uint8', 
                          
                          'dcoilwtico' : 'float64',

                          'const': 'uint8', 
                          'trend': 'uint16',
                          })

test = pd.read_csv('test w features.csv',
                 dtype = {'id': 'uint32', 
                          'store_nbr': 'uint8', 
                          'family': 'category',
                          'onpromotion': 'uint16',

                          'quarter': 'uint8', 
                          'month': 'uint8', 
                          'year': 'uint16',
                          'dayofyear': 'uint16',
                          'dayofmonth': 'uint8',
                          'weekofyear': 'uint8',

                          'Friday' : 'uint8', 
                          'Monday' : 'uint8', 
                          'Saturday' : 'uint8', 
                          'Sunday' : 'uint8', 
                          'Thursday' : 'uint8', 
                          'Tuesday' : 'uint8', 
                          'Wednesday' : 'uint8', 
                          
                          'dcoilwtico' : 'float64',

                          'const': 'uint8', 
                          'trend': 'uint16',
                          })

# oil = pd.read_csv('oil ready.csv')
# holiday_dummies = pd.read_csv('holiday_dummies.csv').fillna(0)

print(train.shape, test.shape)

(3000888, 42) (28512, 41)

test = test.rename(columns={'Unnamed: 0' : 'date'})

print(train.dtypes)
# print(test.dtypes)

date                    object
id                      uint32
store_nbr                uint8
family                category
sales                  float64
onpromotion             uint16
Friday                   uint8
Monday                   uint8
Saturday                 uint8
Sunday                   uint8
Thursday                 uint8
Tuesday                  uint8
Wednesday                uint8
quarter                  uint8
month                    uint8
year                    uint16
dayofyear               uint16
dayofmonth               uint8
weekofyear               uint8
dcoilwtico             float64
const                    uint8
trend                   uint16
sin(1,freq=A-DEC)      float64
cos(1,freq=A-DEC)      float64
sin(2,freq=A-DEC)      float64
cos(2,freq=A-DEC)      float64
sin(3,freq=A-DEC)      float64
cos(3,freq=A-DEC)      float64
sin(4,freq=A-DEC)      float64
cos(4,freq=A-DEC)      float64
sin(5,freq=A-DEC)      float64
cos(5,freq=A-DEC)      float64
sin(6,freq=A-DEC)      float64
cos(6,freq=A-DEC)      float64
sin(7,freq=A-DEC)      float64
cos(7,freq=A-DEC)      float64
sin(8,freq=A-DEC)      float64
cos(8,freq=A-DEC)      float64
sin(9,freq=A-DEC)      float64
cos(9,freq=A-DEC)      float64
sin(10,freq=A-DEC)     float64
cos(10,freq=A-DEC)     float64
dtype: object

def fix_dates(df, datename, set_date=True):
    df['date'] = pd.to_datetime(df[datename], format='%Y-%m-%d')
    df['date'] = df.date.dt.to_period('D')
    if set_date:
        df = df.set_index('date')   
    return df

train = fix_dates(train, 'date')
test = fix_dates(test, 'date')

test.head()

	id	store_nbr	family	onpromotion	Friday	Monday	Saturday	Sunday	Thursday	Tuesday	...	sin(6,freq=A-DEC)	cos(6,freq=A-DEC)	sin(7,freq=A-DEC)	cos(7,freq=A-DEC)	sin(8,freq=A-DEC)	cos(8,freq=A-DEC)	sin(9,freq=A-DEC)	cos(9,freq=A-DEC)	sin(10,freq=A-DEC)	cos(10,freq=A-DEC)
date
2017-08-16	3000888	1	AUTOMOTIVE	0	0	0	0	0	0	0	...	-0.993257	-0.115935	0.796183	-0.605056	-0.154309	0.988023	-0.573772	-0.819015	0.981306	0.192452
2017-08-16	3000889	1	BABY CARE	0	0	0	0	0	0	0	...	-0.993257	-0.115935	0.796183	-0.605056	-0.154309	0.988023	-0.573772	-0.819015	0.981306	0.192452
2017-08-16	3000890	1	BEAUTY	2	0	0	0	0	0	0	...	-0.993257	-0.115935	0.796183	-0.605056	-0.154309	0.988023	-0.573772	-0.819015	0.981306	0.192452
2017-08-16	3000891	1	BEVERAGES	20	0	0	0	0	0	0	...	-0.993257	-0.115935	0.796183	-0.605056	-0.154309	0.988023	-0.573772	-0.819015	0.981306	0.192452
2017-08-16	3000892	1	BOOKS	0	0	0	0	0	0	0	...	-0.993257	-0.115935	0.796183	-0.605056	-0.154309	0.988023	-0.573772	-0.819015	0.981306	0.192452

5 rows × 40 columns

families = test['family'].unique()
stores = test['store_nbr'].unique()

for i in families[:5]:
    print(i)

AUTOMOTIVE
BABY CARE
BEAUTY
BEVERAGES
BOOKS

display(
    train.query('family == @families[1] & store_nbr == @stores[0]').head()
)

	id	store_nbr	family	sales	onpromotion	Friday	Monday	Saturday	Sunday	Thursday	...	sin(6,freq=A-DEC)	cos(6,freq=A-DEC)	sin(7,freq=A-DEC)	cos(7,freq=A-DEC)	sin(8,freq=A-DEC)	cos(8,freq=A-DEC)	sin(9,freq=A-DEC)	cos(9,freq=A-DEC)	sin(10,freq=A-DEC)	cos(10,freq=A-DEC)
date
2013-01-01	1	1	BABY CARE	0.0	0	0	0	0	0	0	...	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000	0.000000	1.000000
2013-01-02	1783	1	BABY CARE	0.0	0	0	0	0	0	0	...	0.103102	0.994671	0.120208	0.992749	0.137279	0.990532	0.154309	0.988023	0.171293	0.985220
2013-01-03	3565	1	BABY CARE	0.0	0	0	0	0	0	1	...	0.205104	0.978740	0.238673	0.971100	0.271958	0.962309	0.304921	0.952378	0.337523	0.941317
2013-01-04	5347	1	BABY CARE	0.0	0	1	0	0	0	0	...	0.304921	0.952378	0.353676	0.935368	0.401488	0.915864	0.448229	0.893919	0.493776	0.869589
2013-01-05	7129	1	BABY CARE	0.0	0	0	0	1	0	0	...	0.401488	0.915864	0.463550	0.886071	0.523416	0.852078	0.580800	0.814046	0.635432	0.772157

5 rows × 41 columns

if zapis:
    for i in train['store_nbr'].unique():
        print(i,'-shape: ', {train.query('store_nbr == @i').shape})
        train.query('store_nbr == @i').to_csv(f'subsets/store_nbr_{i}.csv', index=True) # zapis podzbiorów

Modelowanie

• prognoza będzie polegała na kombinacji dwóch modeli:
  • prostym np. regresja liniowa oparta na celu (sales) i kroku (1 - pierwszy dzień, 1688 - ostatni dzień w przypadku gdy brane byłyby wartości od początku zestawu testowego czyli od 1-01-2013)
  • złożonym np. niezawodny XGBoost oparty na wszystkich pozostałych cechach zestawu
• od faktycznych wartości sprzedaży zostaną usunięte wartości z prostej z regresji liniowej co w wyniku da zdetrendowane wartości sprzedaży, na których będzie uczył się i przewidywał model złożony

Prognozowanie odbywało się osobno dla każdej dziedziny w każdym sklepie, czyli niecałe 1800 modeli, zajmowało to od 12 do 15 minut w zależności od hiperparametrów i złożoności modeli

Sprawdzono kilka sposobów filtrowania dat, ostatecznie zastosowano dwie metody jednocześnie:

• filtr na podstawie eksploracji danych np. ograniczenie zestawu treningowego do wartości tylko po brakach 2015 roku
• filtr na podstawie sumy kumulacyjnej:
  • w przypadkach gdy pierwszy filtr pozostawił daty bez sprzedaży (np. pierwszy filtr mamy od 2015-06-01), a w podzbiorze sprzedaż zaczyna się od 2016-01-01 i mamy kilka miesięcy braków
  • przy pomocy sumy kumulacyjnej można sprawdzić przy której dacie zaczyna się rzeczywista sprzedaż - dla powyższego przykładu od 2015-06-01 do 2015-12-31 suma kumulacyjna wynosiłaby 0, a w 2016-01-01 już np. 40
  • kolejno przy pomocy maski wartości True/False czy wartość jest nierówna 0 można odfiltrować puste daty

Wczytanie danych z optymalnym formatem:

test = pd.read_csv('test w features.csv',
                 dtype = {'id': 'uint32', 
                          'store_nbr': 'uint8', 
                          'family': 'category',
                          'onpromotion': 'uint16',

                          'quarter': 'uint8', 
                          'month': 'uint8', 
                          'year': 'uint16',
                          'dayofyear': 'uint16',
                          'dayofmonth': 'uint8',
                          'weekofyear': 'uint8',

                          'Friday' : 'uint8', 
                          'Monday' : 'uint8', 
                          'Saturday' : 'uint8', 
                          'Sunday' : 'uint8', 
                          'Thursday' : 'uint8', 
                          'Tuesday' : 'uint8', 
                          'Wednesday' : 'uint8', 
                          
                          'dcoilwtico' : 'float64',

                          'const': 'uint8', 
                          'trend': 'uint16',
                          })

test['date'] = pd.to_datetime(test['Unnamed: 0'], format='%Y-%m-%d')
test['date'] = test.date.dt.to_period('D')
test = test.set_index('date')
test = test.drop('Unnamed: 0', axis=1)

print(test.shape)
test.head(3)

(28512, 40)

	id	store_nbr	family	onpromotion	Friday	Monday	Saturday	Sunday	Thursday	Tuesday	...	sin(6,freq=A-DEC)	cos(6,freq=A-DEC)	sin(7,freq=A-DEC)	cos(7,freq=A-DEC)	sin(8,freq=A-DEC)	cos(8,freq=A-DEC)	sin(9,freq=A-DEC)	cos(9,freq=A-DEC)	sin(10,freq=A-DEC)	cos(10,freq=A-DEC)
date
2017-08-16	3000888	1	AUTOMOTIVE	0	0	0	0	0	0	0	...	-0.993257	-0.115935	0.796183	-0.605056	-0.154309	0.988023	-0.573772	-0.819015	0.981306	0.192452
2017-08-16	3000889	1	BABY CARE	0	0	0	0	0	0	0	...	-0.993257	-0.115935	0.796183	-0.605056	-0.154309	0.988023	-0.573772	-0.819015	0.981306	0.192452
2017-08-16	3000890	1	BEAUTY	2	0	0	0	0	0	0	...	-0.993257	-0.115935	0.796183	-0.605056	-0.154309	0.988023	-0.573772	-0.819015	0.981306	0.192452

3 rows × 40 columns

stores = test['store_nbr'].unique()
families = test['family'].unique()

families

['AUTOMOTIVE', 'BABY CARE', 'BEAUTY', 'BEVERAGES', 'BOOKS', ..., 'POULTRY', 'PREPARED FOODS', 'PRODUCE', 'SCHOOL AND OFFICE SUPPLIES', 'SEAFOOD']
Length: 33
Categories (33, object): ['AUTOMOTIVE', 'BABY CARE', 'BEAUTY', 'BEVERAGES', ..., 'PREPARED FOODS', 'PRODUCE', 'SCHOOL AND OFFICE SUPPLIES', 'SEAFOOD']

holidays = pd.read_csv('holiday_dummies.csv')
holidays['date'] = pd.to_datetime(holidays['date'], format='%Y-%m-%d')
holidays['date'] = holidays.date.dt.to_period('D')
holidays = holidays.set_index('date')
holidays = holidays.astype('uint8')
holidays.dtypes

description_Batalla de Pichincha                      uint8
description_Black Friday                              uint8
description_Carnaval                                  uint8
description_Cyber Monday                              uint8
description_Dia de Difuntos                           uint8
                                                      ...  
description_Traslado Batalla de Pichincha             uint8
description_Traslado Independencia de Guayaquil       uint8
description_Traslado Primer Grito de Independencia    uint8
description_Traslado Primer dia del ano               uint8
description_Viernes Santo                             uint8
Length: 72, dtype: object

Funkcja wczytywanie podzbiorów treningowych

def get_subset(store_number, family, filter,): # family jest w .query()
    current_store = pd.read_csv(f'subsets\store_nbr_{store_number}.csv', # odpowiedni format
                dtype = {'id': 'uint32', 
                        'store_nbr': 'uint8', 
                        'family': 'category',
                        'onpromotion': 'uint16',

                        'quarter': 'uint8', 
                        'month': 'uint8', 
                        'year': 'uint16',
                        'dayofyear': 'uint16',
                        'dayofmonth': 'uint8',
                        'weekofyear': 'uint8',

                        'Friday' : 'uint8', 
                        'Monday' : 'uint8', 
                        'Saturday' : 'uint8', 
                        'Sunday' : 'uint8', 
                        'Thursday' : 'uint8', 
                        'Tuesday' : 'uint8', 
                        'Wednesday' : 'uint8', 
                        
                        'dcoilwtico' : 'float32',

                        'const': 'uint8', 
                        'trend': 'uint16',
                        })
    
    current_store = current_store.query('family == @family').copy()

    current_store['date'] = pd.to_datetime(current_store['date'], format='%Y-%m-%d') # poprawka dat
    current_store['date'] = current_store.date.dt.to_period('D')
    current_store = current_store.set_index('date')

    current_store = current_store.join(holidays, how='left').fillna(0) # złączenie świąt

    # filtr sumy kumulacyjnej, z poprawką na dziedziny bez żadnej sprzedaży
    if current_store['sales'].sum() != 0:
        mask = current_store['sales'].cumsum().ne(0)
        current_store = current_store[mask]



    return current_store.query('index >= @filter') # filtr po datach z eksploracji danych

for i in stores[[0]]:
    for j in families[[1]]:
        current_subset = get_subset(i, j, '2013-09-01')

current_subset

	id	store_nbr	family	sales	onpromotion	Friday	Monday	Saturday	Sunday	Thursday	...	description_Terremoto Manabi+5	description_Terremoto Manabi+6	description_Terremoto Manabi+7	description_Terremoto Manabi+8	description_Terremoto Manabi+9	description_Traslado Batalla de Pichincha	description_Traslado Independencia de Guayaquil	description_Traslado Primer Grito de Independencia	description_Traslado Primer dia del ano	description_Viernes Santo
date
2013-09-01	433027	1	BABY CARE	0.0	0	0	0	0	1	0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
2013-09-02	434809	1	BABY CARE	0.0	0	0	1	0	0	0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
2013-09-03	436591	1	BABY CARE	0.0	0	0	0	0	0	0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
2013-09-04	438373	1	BABY CARE	0.0	0	0	0	0	0	0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
2013-09-05	440155	1	BABY CARE	0.0	0	0	0	0	0	1	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
2017-08-11	2991979	1	BABY CARE	0.0	0	1	0	0	0	0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0
2017-08-12	2993761	1	BABY CARE	0.0	0	0	0	1	0	0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
2017-08-13	2995543	1	BABY CARE	0.0	0	0	0	0	1	0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
2017-08-14	2997325	1	BABY CARE	0.0	0	0	1	0	0	0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
2017-08-15	2999107	1	BABY CARE	0.0	0	0	0	0	0	0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

1445 rows × 113 columns

Funkja modeli

from sklearn.linear_model import LinearRegression
from xgboost import XGBRegressor

def hybrid_model(train, test, plot=False):
    
    ## LINEAR REGRESSION - TREND
    model_lr = LinearRegression()
    model_lr.fit(train[['trend']], train[['sales']])

    lr_fit_pred = model_lr.predict(train[['trend']])
    lr_fit_pred = pd.DataFrame(lr_fit_pred, index=train.index)


    ## ręczna kontynuacja prostej regresji liniowej, coś się krzaczyło
    step = lr_fit_pred[0][-1] - lr_fit_pred[0][-2]
    start = lr_fit_pred[0][-1] + step
    amount = 16

    manual_lr = np.linspace(start, start+(step*amount),amount,endpoint=False) 
    lr_fore_pred = manual_lr
    lr_fore_pred = pd.DataFrame(lr_fore_pred, index=test.index)

    ## detrend
    train_detrend = train.copy()
    train_detrend['sales'] =  train_detrend['sales'] - lr_fit_pred[0]


    ## XGBoostRegressor - reszta cech na zdetrendowanych danych
    model_xgb = XGBRegressor(n_estimators=25)
    X_detrend_train = train_detrend.iloc[:, 4:]
    y_detrend_train = train_detrend['sales']

    model_xgb.fit(X_detrend_train, y_detrend_train)
    detrend_test_fore = model_xgb.predict(test.iloc[:, 3:])

    test_fore = pd.DataFrame(detrend_test_fore + lr_fore_pred[0], index=test.index,).clip(0.0)
    test_fore = test_fore.rename(columns={0:'Sales'})
    test_fore['Id'] = test['id']

    ## sprawdzenie czy jest sens zaokrąglić np. sprzedaż książek to będą tylko liczby naturalne, jeśli tak - zaokrąglamy
    if train['sales'].apply(float.is_integer).all():
        test_fore["Sales"] = test_fore["Sales"].round(0)


    finished_result = test_fore.reset_index().set_index('Id').drop('date', axis=1)        

    return [lr_fit_pred, lr_fore_pred, train_detrend, test_fore, finished_result.astype('float16')]

Filtr dat

date_filter = {
'AUTOMOTIVE' :                  '2013-01-01', # normalne 
'BABY CARE' :                   '2015-12-07', 
'BEAUTY' :                      '2013-01-01', # normalne 
'BEVERAGES' :                   '2015-05-28', 
'BOOKS' :                       '2017-01-01', # normalne 
'BREAD/BAKERY' :                '2013-01-01', # normalne 
'CELEBRATION' :                 '2015-06-01', 
'CLEANING' :                    '2013-01-01', # normalne 
'DAIRY' :                       '2013-11-07', 
'DELI' :                        '2013-01-01', # normalne 
'EGGS' :                        '2013-01-01', # normalne 
'FROZEN FOODS' :                '2013-01-01', # normalne 
'GROCERY I' :                   '2013-01-01', # normalne 
'GROCERY II' :                  '2013-01-01', # normalne 
'HARDWARE' :                    '2013-01-01', # normalne 
'HOME AND KITCHEN I' :          '2014-09-01', 
'HOME AND KITCHEN II' :         '2015-06-01', 
'HOME APPLIANCES' :             '2013-01-01', # normalne 
'HOME CARE' :                   '2015-05-04', 
'LADIESWEAR' :                  '2015-06-01', 
'LAWN AND GARDEN' :             '2016-12-03', 
'LINGERIE' :                    '2013-01-01', # normalne 
'LIQUOR,WINE,BEER' :            '2016-05-10', 
'MAGAZINES' :                   '2015-09-30', 
'MEATS' :                       '2016-10-07', 
'PERSONAL CARE' :               '2013-01-01', # normalne 
'PET SUPPLIES' :                '2015-06-01', 
'PLAYERS AND ELECTRONICS' :     '2015-06-01', 
'POULTRY' :                     '2013-11-04', 
'PREPARED FOODS' :              '2013-01-01', # normalne 
'PRODUCE' :                     '2015-06-01', 
'SCHOOL AND OFFICE SUPPLIES' :  '2013-01-01', # normalne 
'SEAFOOD' :                     '2013-01-01', # normalne 
}

Pętla prognozowania

import random
plot_q = False # czy wykreślać i zapisywać losową część prognoz

# podstawa ramki danych do dodawania kolejnych prognoz, na koniec zostaną posortowane po `ID` przed wrzuceniem wyników na Kaggle
base = pd.DataFrame(columns=['Sales', 'Id'])
base = base.set_index('Id') 

for i in stores[[0]]: # na potrzeby tego zeszytu w portfolio tylko jeden sklep i jedna dziedzina zostanie zaprognozowana
    for j in families[[20]]:
        filter = date_filter[j] # która data 
        current_subset = get_subset(i, j, filter)

        current_test = test.query('store_nbr == @i & family == @j')
        current_test = current_test.join(holidays, how='left').fillna(0)

        result = hybrid_model(current_subset, current_test)

        base = pd.concat([base, result[4].astype('float16')])        

        # if random.randint(1, 100) == 1:
        # if plot_q:
        if True:
            name = j.capitalize().replace("/", " ")

            fig = go.Figure()
            fig.add_trace(go.Scatter(x=current_subset.index.to_timestamp()[-365:], y=current_subset['sales'][-365:], name=f'{j} sales')) 
            fig.add_trace(go.Scatter(x=result[3].index.to_timestamp(), y=result[3]['Sales'], mode='lines', name=f'{j} forecast'))
            fig.add_trace(go.Scatter(x=result[0].index.to_timestamp()[-365:], y=result[0][0][-365:], mode='lines', name=f'LR fit'))
            fig.add_trace(go.Scatter(x=result[1].index.to_timestamp(), y=result[1][0], mode='lines', name=f'LR fore'))

            fig.update_layout(title=f'Forecast for {name} - Store {i}', 
                            width=1200, height = 400,
                            template='plotly_dark',
                            legend=dict(
                                orientation="h",
                                yanchor="bottom",
                                y=1.02,
                                xanchor="right",
                                x=1))            
            pio.write_image(fig, f"subset images 4/{i} {name}.png")
            
            pio.write_image(fig, f"subset images cover/{i} {name}.png")



check_png(fig)

base = base.reset_index().drop_duplicates(subset='Id').set_index('Id').sort_index()
base

	Sales
Id
3000908	16.0
3002690	9.0
3004472	13.0
3006254	17.0
3008036	17.0
3009818	18.0
3011600	16.0
3013382	24.0
3015164	19.0
3016946	19.0
3018728	14.0
3020510	17.0
3022292	16.0
3024074	17.0
3025856	25.0
3027638	19.0

# sprawdzenie czy rzeczywiście nie ma powtórzeń
print(base.shape)
base.index.value_counts().max()

(16, 1)

1

# tu mamy ograniczoną pętlę do tylko jednego sklepu i jednej dziedziny, w pełnej petli mamy odpowiednią ilość
print(16 * 54 * 33)
print(test.shape)

28512
(28512, 40)

if zapis:
    base.to_csv('my submission 4 - LR, 12 est, manual and cumsum date filter.csv')

Wnioski

Ogółem rzecz biorąc, prognozowanie wyżej przedstawioną metodą okazało się dość trafne - spośród wszystkich podejść hybryda regresji liniowej i dość prostego XGBoost'a (12 estymatorów) dała wyniki RMSLE (Root Mean Squared Logarithmic Error) równy 0.46, zajmując 88 miejsce spośród 681 uczestników (13 percentyl). W gwoli ścisłości RMSLE liczony jest wg wzoru:

$$\sqrt{ \frac{1}{n} \sum_{i=1}^n \left(\log (1 + \hat{y}_i) - \log (1 + y_i)\right)^2}$$

gdzie:

• $\hat y_i$ to przewidywana wartość
• $y_i$ to rzeczywista wartość
• $n$ to ilość obserwacji

Poniżej przedstawiono tabelę z wynikami różnych podejść:

nr podejścia	Wynik	Pozycja	Opis
1	0.50463	122	LR + XGB, manual date filter
2	0.49233	114	LR + XGB, cumsum() date filter
3	0.50251		LR + XGB (333 estimators), cumsum() date filter
4	0.48187	108	LR + XGB (50 estimators), cumsum() date filter
5	0.47315	98	LR + XGB (25 estimators), cumsum() date filter
6	0.46249	88	LR + XGB (12 estimators), manual and cumsum() date filter
6	0.46807		LR + XGB (5 estimators), manual and cumsum() date filter
7	0.47727		Ridge $\alpha=1$ + XGB (25 estimators), manual and cumsum() date filter
8	0.46909		Ridge $\alpha=3$ + XGB (12 estimators), manual and cumsum() date filter
9	0.46909		Ridge $\alpha=9$ + XGB (12 estimators), manual and cumsum() date filter

Okazało się, że zaawansowany model w hybrydzie powinien być dość prosty - najlepsze wyniki przy 12 drzewach (estymatorach) w XGBoost. Regresja grzbietowa (Ridge Regression) zmniejszająca dopasowanie do danych treningowych nie przyniosła lepszych efektów.

Zalecenia jak usprawnić model:

• dokładniejsza eksploracja danych - to jest niestety trudne w przypadku podejścia, gdzie każdy sklep i jego dziedziny sprzedaży modeluje się osobno
• lepsze przyjrzenie się świętom - prawdopodobnie wrzucenie wszystkich świąt niepotrzebnie powoduje szumy, na podstawie eksploracji można by zaznaczyć kilka świąt / wydarzeń - Wielkanoc, początek roku szkolnego (sprzedaż artykułów biurowych i szkolnych), Święta Bożego Narodzenia, Nowy Rok, istotniejsze święta narodowe Ekwadoru
• zastosowanie innych modeli np. regresji LASSO, Elastic Net, sieci neuronowych
• zastosowanie zupełnie innego podejścia - np. modelowania ARIMA specjalnie przeznaczonego do serii czasowych
• zastosowanie rekursywnej metody predykcji:
  • każda pojedyncza predykcja (sprzedaż danego dnia, danej dziedziny w danym sklepie) zostaje dołączona do zestawu treningowego i kolejno następuje ponowne nauczanie modelu i kolejna pojedyncza predykcja
  • zaletą tej metody mogą być celniejsze wyniki, z drugiej strony wczesny błąd w predykcjach będzie "niesiony" z każdą kolejną predykcją,
  • metodę tę można poprawić dodając lagi wartości, którą chcemy przewidzieć (tu już pojawiają się zagadnienia z ekonometrii), lagi to wartości z przeszłości
  • główną wadą metody rekursywnej są wysokie wymagania obliczeniowe